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Abstract. DMRT-ML is a physically based numerical model 
designed to compute the thermal microwave emission of 
a given snowpack. Its main application is the simula- 
tion of brightness temperatures at frequencies in the range 
1-200 GHz similar to those acquired routinely by space- 
based microwave radiometers. The model is based on the 
Dense Media Radiative Transfer (DMRT) theory for the 
computation of the snow scattering and extinction coeffi- 
cients and on the Discrete Ordinate Method (DISORT) to nu- 
merically solve the radiative transfer equation. The snowpack 
is modeled as a stack of multiple horizontal snow layers and 
an optional underlying interface representing the soil or the 
bottom ice. The model handles both dry and wet snow con- 
ditions. Such a general design allows the model to account 
for a wide range of snow conditions. Hitherto, the model 
has been used to simulate the thermal emission of the deep 
fim on ice sheets, shallow snowpacks overlying soil in Arc- 
tic and Alpine regions, and overlying ice on the large ice- 
sheet margins and glaciers. DMRT-ML has thus been vali- 
dated in three very different conditions: Antarctica, Barnes 
Ice Cap (Canada) and Canadian tundra. It has been recently 
used in conjunction with inverse methods to retrieve snow 
grain size from remote sensing data. The model is written in 
Fortran90 and available to the snow remote sensing commu- 
nity as an open-source software. A convenient user interface 
is provided in Python. 


1 Introduction 

Passive microwave radiometers on-board satellites acquire 
useful observations for the characterization of snow-covered 
areas. These observations are available in these areas several 
times a day, are relatively insensitive to the atmosphere in 
many frequency bands, and are independent of the solar il- 
lumination. They are sensitive to several properties relevant 
for monitoring the snow cover and have been exploited in nu- 
merous algorithms to retrieve continental snow cover extent 
(Grody and Basist, 1996), snow depth and snow water equiv- 
alent on both land (Josberger and Mognard, 2002; Kelly and 
Chang, 2003; Derksen et ah, 2003) and sea ice (Cavalieri 
et ah, 2012; Brucker and Markus, 2013), snow accumulation 
on ice sheets (Abdalati and Steffen, 1998; Vaughan et ah, 
1999; Winebrenner et ah, 2001; Arthern et ah, 2006), melt 
events (Abdalati and Steffen, 1997; Picard and Fily, 2006), 
snow temperature (Shuman et ah, 1995; Schneider, 2002; 
Schneider et ah, 2004), and snow grain size (Brucker et ah, 
2010; Picard et ah, 2012). Some of these studies are based 
on empirical relationships supported by physical interpre- 
tations (Koenig et ah, 2007) and others directly use phys- 
ical models and data assimilation techniques (Durand and 
Margulis, 2007; Picard et ah, 2009; Takala et ah, 201 1; Toure 
et ah, 201 1; Huang et ah, 2012). In both cases, understanding 
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and modeling the physical processes of the microwave emis- 
sion by snow and the underlying surface are crucial. 

Modeling the snow microwave emission from snow phys- 
ical properties takes in general two successive steps. The 
first step is the computation of the electromagnetic proper- 
ties (e.g., effective dielectric constant, scattering and absorp- 
tion coefficients) that characterize the propagation and the 
single events of interaction between wave and matter. These 
properties are calculated from the geometric micro-structural 
properties of the snow, assumed to be homogeneous within a 
given layer. The second step is the computation of the emis- 
sion and the propagation through the snowpack by account- 
ing for the multiple interactions of microwaves within the 
snow as well as by refraction, reflection and transmission at 
the interfaces. It is often treated with the radiative transfer 
equation on a plane-parallel medium for which generic solu- 
tions exist (Tsang et al., 1985; Fung, 1994). In contrast, the 
first step is specific to the medium. In the case of snow, the 
main challenge for the electromagnetic calculation is the high 
density of scatterers. The volume occupied by the scatterers 
over the total volume (referred to as the fractional volume /) 
is significant - typically larger than 20 % which implies 
that the scatterers strongly interact with each other and the 
independent scattering theory used for vegetation or clouds 
(Tsang et al., 1985; Ulaby et al., 1986; Chuah and Tan, 1989) 
is inadequate. Several empirical formulas relating the scat- 
tering and absorption coefficients to grain size and density 
have been proposed to solve this issue and are found in the 
Helsinki University of Technology model (HUT/TKX) (Pul- 
liainen et al., 1999; Lemmetyinen et al., 2010) and in the Mi- 
crowave Emission Model of Layered Snowpacks (MEMLS) 
(Wiesmann and Matzler, 1999). The derivation of relation- 
ships from Maxwell’s equations is an attractive alternative 
because it is independent of particular snow conditions and 
the assumptions are explicit: 

- the strong fluctuation theory (SFT, Stogryn, 1986) al- 
lows calculations at low frequencies typically less than 
20 GHz; 

- to cover the frequency range of 1-100 GHz, which is 
relevant for existing radiometers, Matzler (1998) pro- 
posed the improved Born approximation (IBA). As in 
the SFT, the size of grains is given by the correlation 
length (Matzler, 2002). This metric is well defined but 
direct measurements are only possible from 2-D or 3- 
D micro-structure data, which require advanced experi- 
mental techniques (Wiesmann et al., 1998). Estimations 
can be indirectly obtained from quantities measurable 
in the field like the snow specific surface area (Matzler, 
2002; Arnaud et al., 201 1) or the micro-penetration pro- 
file (Lowe and van Herwijnen, 2012); 

- the dense media radiative transfer theory (DMRT; 
Tsang et al., 1985; West et al., 1993; Shih et al., 1997; 
Tsang et al., 2000a, 2007) considers snow as a collec- 


tion of spherical ice particles and provides analytical 
formulas of the effective propagation wave number, the 
scattering and absorption coefficients as a function of 
the sphere radius and the density. It is attractive be- 
cause sphere radius and snow density are well defined 
and, when the snow grains are close to spheres like in 
fine grained snow (Colbeck, 1993), they are both mea- 
surable in the field. Moreover, the DMRT-ML theory 
has solid theoretical grounds and has been regularly im- 
proved during the last two decades. 

These theories have been implemented in several models: 
the SFT in Surdyk and Fily (1995), the IBA in MEMLS 
(Matzler and Wiesmann, 1999) and the DMRT-ML in Ma- 
celloni et al. (2001); Tedesco and Kim (2006); Liang et al. 
(2008); Grody (2008); Brogioni et al. (2009). In addition to 
the underlying theory, the models differ by many aspects, in- 
cluding the methods for solving the radiative transfer, the 
presence of smooth or rough interfaces, the possible num- 
ber of layers, etc. In particular, HUT uses the two-stream 
method and MEMLS the six-stream method while most 
DMRT-based models consider a larger number of streams. 
The comparison for a large variety of snow types (Tedesco 
and Kim, 2006) showed that no particular model systema- 
tically reproduces all of the experimental data. It is not yet 
known if the discrepancies were attributable to the electro- 
magnetic theory, specific details of the implementations of 
the models, or uncertainties pertaining to input or evalua- 
tion data. Moreover, the different representations of the snow 
grain size in these different approaches are a major limit to 
the comparison. Both MEMLS and HUT/TKK models are 
widely used (Butt and Kelly, 2008; Durand et al., 2008; Rees 
et al., 2010; Brucker et al., 2011b; Gunn et al., 2011). In con- 
trast, several groups use home-made DMRT-based models 
but no widely-used reference implementation exists, which 
limits the spread of this theory and limits the comparisons 
between studies. 

The objective of this paper is to present the DMRT-ML 
(DMRT Multi-Layer) model that was released under an open 
source license and accompanied by a detailed documentation 
(http://lgge.osug.fr/~picard/dmrtml/). This model initially 
developed at Laboratoire de Glaciologie and Geophysique de 
(’Environment in Grenoble, France, for perennial snowpack 
(Brucker et al., 2010, 2011a) and improved for seasonal snow 
on soil (Roy et al., 2013) and on superimposed ice (Dupont 
et al., 2012), is now suitable for modeling the microwave 
emission in a wide range of snowy continental environments. 
It is also designed to be extensible and uses standard For- 
tran90, which allows efficient computations with different 
types of computers and operating systems and facilitates the 
embedding in other models or assimilation schemes. A con- 
venient user interface is provided by the optional bindings 
with the Python language. 
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This paper is a comprehensive scientific reference for 
DMRT-ML users. It is structured as follows. Section 2 
describes DMRT-ML in detail. Section 3 presents the 
sensitivity to the most important input variables and parame- 
ters and provides practical recommendations on the validity 
range of the input variables and parameters to use. Section 4 
summarizes the results of the detailed comparisons (Brucker 
et al., 2011a; Roy et al., 2013; Dupont et al., 2012) with ra- 
diometric data obtained for various cold environments. 


2 The DMRT-ML model 

2.1 Radiative transfer equation and model architecture 

The energy emanating from snowpacks is the result of the 
thermal emission of the snow, the substratum and the atmo- 
sphere, and of the complex propagation of this energy toward 
the upper snow layers. The emission and propagation can be 
described by the radiative transfer equation. Assuming that 
the medium is a stack of L plane-parallel layers containing 
an isotropic and homogeneous material, the equation in every 
layer is (Jin, 1994) 

cos0-^-Tb (z,0,4>) = -/c e T B (z,0,(t>) 

d z 

nil In 

+ J J sm9'd9'd(P'V(9,(l),9',4>')TB(z,9\(P') + K :i T(z)I. (l) 

0 0 

By convention, matrices are written in bold, and vectors 
in bold and italic. I is the unit column vector. T B (z,0,cp) 
is the radiance at depth z propagating along the direction 
with zenith angle 6 and azimuth angle <fi. According to the 
Rayleigh Jeans approximation, that is valid in the microwave 
range, T B (z,6,<f>) also represents the brightness tempera- 
ture. Since the medium is isotropic, T B is a vector with only 
two non-null components (Jin, 1994, p. 19): the vertically 
and horizontally polarized brightness temperature. K e and /c a 
are the extinction and absorption coefficients in the layer and 
P is the phase function. T is the physical temperature of the 
layer. 

The model solves this equation for the particular medium 
described by the input variables and parameters (Fig. 1). 
Several steps are needed: the determination of the dielectric 
properties of the constitutive materials (Sect. 2.2), the de- 
termination of the effective propagation constant and of the 
extinction and scattering coefficients for each homogeneous 
layer of the medium (Sect. 2.3) by using the DMRT-ML the- 
ory, the determination of the boundary conditions (Sect. 2.4) 
and lastly the numerical computation of the solution of the ra- 
diative transfer equation (Sect. 2.5). The result is the bright- 
ness temperature emerging from the surface in several direc- 
tions. Finally, other relevant variables are accessible via ad- 
ditional post-computations (Sect. 2.6). 
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Fig. 1. The snowpack viewed by DMRT-ML. L is the number of 
layers. For each homogeneous layer k, the input variables and pa- 
rameters are grain size a k , snow density p k , physical temperature 
T k , stickiness r k and liquid water content LWC&. The extinction 
and absorption coefficients K t and k & are calculated with the DMRT- 
ML theory. The boundary conditions require the input of the down- 
welling brightness temperature latinos and the variables and param- 
eters to compute the reflection coefficient r(9) of the underlying 
interface (see Table 1). The DISORT method is used to solve the ra- 
diative transfer equation in the snowpack and provides as output the 
brightness temperature T B (9, z = 0) emerging from the snowpack. 


2.2 Dielectric constants 


The computation of the scattering and extinction coefficients 
in the DMRT-ML theory requires the dielectric constant of 
the scatterer and background materials. When snow is dry 
(temperature strictly below the melting point), the grains are 
assumed to be composed of pure ice, whose dielectric con- 
stant £; ce was measured by Matzler and Wegmiiller (1987) 
and Matzler et al. (2006): 


e ice = 3.1884 + 0.00091(7 - 273.0) + j + /3v) . 
a = (0.00504 + 0.00620) exp(— 22.1©), 

0.0207 


exp(^) 

(exp(^)-iy 


(2a) 

(2b) 

(2c) 


+ 1.1610 -11 v 2 + exp (-9.963 +0.0372(7 - 273.16)) , 

(2d) 

300 

0= 1, (2e) 


where v is the frequency in GHz, 7 the ice temperature 
in K and j 2 — —1. This model is valid for temperatures 
above 240 K and in the microwave range (1-200 GHz) but 
its accuracy certainly decreases at low frequencies (Jiang and 
Wu, 2004). 
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When the snow is wet (at 273.15 K), ice grains are coated 
by liquid water with an uneven distribution controlled by 
capillarity forces. Since such a complex particle cannot be 
taken explicitly into account in the DMRT-ML theory and 
given the low content of liquid water usually present in moist 
and wet snow (less than 8%; Colbeck, 1993), DMRT-ML 
assumes that the grains are composed of a fictitious homo- 
geneous material. The dielectric constant £ wet ; ce of this ma- 
terial is calculated using the mixture relationship (Borghese 
et ah, 2007; Chopra and Reddy, 1986): 


f-'wct ice — 


C+ + 2 C- 
C+-C- 


£ water , 


C+ — £ice T" £water 


C— — (fiice £ water)(l LWC), 


(3a) 

(3b) 

(3c) 


where LWC is the liquid water content expressed here as the 
ratio between the volume of liquid water and the volume of 
ice present in snow (in m 3 m -3 ) and the water dielectric con- 
stant £ w ater is given by 


, et — £2 , £0 — £l 

£water — £2 + “ y r “ p" , 
J V 2 ^ U | 

£ 0 = 77.66 + 103.30, 

£ i = 0.067l£o , 

£2 = 3.52 + 7.520, 

vi — 20.2 — 146.40 -3160 2 , 

r>2 = 39. 8vj . 


(4a) 

(4b) 

(4c) 

(4d) 

(4e) 

(4f) 


2.3 Snow extinction and absorption coefficients, and 
phase function: the DMRT-ML theory 


The extinction and absorption coefficients, as well as the 
form of the phase function are obtained by the so-called 
DMRT-ML theory. Different versions of this theory have 
been published over time and the underlying assumptions 
can significantly differ from each other. Four characteristics 
distinguish these versions: (i) the underlying approximation 
used for the DMRT-ML derivation, (ii) the particle size with 
respect to the wavelength, (iii) the stickiness between parti- 
cles and (iv) the distribution of particle size. DMRT-ML pro- 
poses some of these versions and allows various options as 
detailed here. 

DMRT-ML uses the quasicrystalline approximation with 
coherent potential (QCA-CP) (Tsang et ah, 2000b) and is 
limited to particle size smaller than the wavelength. This 
may be a limitation at frequencies higher than 37 GFIz and 
with large grains commonly found in aged snow. The cal- 
culation for large particles requires a Mie-like development 
(Tsang et ah, 2000a) and is computationally much more in- 
tensive than the QCA-CP calculation. In addition, it leads to a 
form of the phase function that is incompatible with the opti- 
mization of the radiative transfer equation resolution used in 
DMRT-ML (Sect. 2.5). To avoid the strong divergence that 


characterized QCA-CP with large particles, Grody (2008) 
proposed an empirical and computationally efficient treat- 
ment of this issue. He noticed that snow is composed of par- 
ticles with a broad range of sizes, which results in a smooth- 
ing of the undulation characteristic of the Mie resonances. 
Hence, a good approximation of the medium scattering ef- 
ficiency <2s f°r large particles is the asymptotic limit, i.e., 
<2 S ^ 2, accounting from the fact that the absorption is weak 
(Twomey and Bohren, 1980). If enabled by the user, DMRT- 
ML applies this idea by limiting the Q s maximum value to 2. 
Nevertheless, this ad hoc correction has not the objective to 
replace the rigorous solution for fine-grained studies (Tsang 
et ah, 2000a). 

Recent versions of DMRT-ML introduce the concept of 
stickiness. Instead of considering randomly positioned non- 
penetrable spheres, the sticky spheres are attracted to one an- 
other and tend to form clusters with large voids between. This 
concept is meant to better represent the micro-structure of 
natural snow using solely spherical grains. In this case, the 
stickiness is also a means to account for coarse grained snow 
by considering that such snow is made of small clustered 
particles. However, DMRT-ML only implements the “short 
range” version of the sticky formulation, which assumes that 
the clusters are small with respect to the wavelength (Tsang 
et ah, 2000b, pp. 504-505). In this simplified case, the phase 
function remains identical to that of the nonsticky small par- 
ticle case for which the optimization of the resolution of the 
radiative transfer equation works (Sect. 2.5). 

Another way to improve the representation of the micro- 
structure of natural snow is to consider a collection of spheres 
with different sizes. In DMRT-ML, the particle sizes follow 
the Rayleigh distribution (Jin, 1994). In addition, only the 
nonsticky particle case is available because the formulation 
with both stickiness and size distribution leads to a quadratic 
system of equations that is difficult to solve and is computa- 
tionally intensive (Tsang and Kong, 2001, p. 430). 

In summary, DMRT-ML includes two implementations: 

- QCA-CP mono-disperse, with optional “short range” 
stickiness, and optional Grody’s-based empirical treat- 
ment of large particles; 

- QCA-CP poly-disperse with a Rayleigh distribution, no 
stickiness and no large particles. 

The mono-disperse version is formulated according to 
Shih et al. (1997). The effective dielectric constant with- 
out scattering £ e fl 0 is obtained by solving the follow- 
ing quadratic Eq. (3) in Shih et al. (1997) with a — 0 or 
Eq. (5.3.125) in Tsang and Kong (2001): 

Eli TO + £efffl - 4/) - e b ) - Cb^-d -/) = «, (5) 

where / is the fractional volume of scatterers, e b and e s are 
the dielectric constants of the background and scatterers, re- 
spectively. The effective dielectric constant with scattering 
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(Eq. 3 in Shih et al., 1997) combined with Eq. (5) yields: 


E e & — £b + (£effl) - £b) 

+ 7^( A:ofl )V £ ' e fiO £s t £b 
, gs-£b „ (l-/) 4 \ 

+ 3£ e ffl) /) (l+2/-t/(l-/)) 2 j’ 

where a is the radius of the spheres and ko — jt jX is the wave 
number with X the wavelength, t is zero for nonsticky spheres 
and otherwise is the largest solution of Eq. 6 in Shih et al. 
(1997) 


/ 

12 


(t + 


/ 

1 -/ 


)t + 


1 + // 2 
(1 - f) 2 


— 0, 


(7) 


where r is the stickiness parameter (Shih et al., 1997; Tsang 
and Kong, 2001, p. 430). At last, the extinction and absorp- 
tion coefficients are respectively given by 


Ke = 2A: 0 3v / £ ; eff, 

( l -/) 4 

(l+2/-r/(l-/)) 2 ’ 


*s = yko a 3 f 




i+Sd-/) 


(8a) 

(8b) 

(8c) 


with 3 the imaginary part indicator and k s the scattering co- 
efficient. 

In the case of small particles, the phase function in 
DMRT-ML has the same form as with independent particles 
Rayleigh phase function (Jin, 1994, Eqs. 2-26 to 2-28). 

The DMRT-ML theory inherits its name from the fact that 
it extends the classical radiative transfer theory, which re- 
quires the particles not to interact with each other. This con- 
dition is met only when the fractional volume / is less than 
a few percent (Ishimaru and Kuga, 1982). Despite its name, 
recent works suggest the DMRT-ML theory does not work 
for any large fractional volumes (see Sect. 3.3 for details). 
As a consequence, for icy layers or layers subject to several 
melt-refreeze cycles, it is recommended to use the DMRT- 
ML theory for a collection of air bubbles embedded in ice 
background. This is achieved in DMRT-ML by exchanging 
the dielectric constant of ice and air in Eqs. (5) and (6). 


2.4 Properties of the interfaces 

Once the extinction and absorption coefficients and the phase 
function within every layer are determined, all the variables 
and parameters in Eq. (1) are known and a general solu- 
tion can be obtained independently for every layer. To obtain 
the particular solution and hence the brightness temperature 
emerging from the surface, it is necessary to propagate the 
radiation between the individual snow layers. This propaga- 
tion must ensure the energy conservation and requires the 
reflection properties at every interface as well as any external 
source of energy. 


For the interfaces within the snowpack and at the air- 
snow interface, DMRT-ML considers smooth interfaces. The 
roughness of this interface due to ripples, sastrugi and other 
dunes (Watanabe, 1978) may contribute especially at hori- 
zontal polarization and grazing incidence angles (Lacroix 
et al., 2009), but the existing formulations of the bi-static re- 
flection coefficient are based on the method of moments and 
are therefore computationally intensive and limited to low 
frequencies (< 18 GHz; Liang et al., 2009; Chang and Tsang, 
2011). For flat interfaces, the reflection coefficient only de- 
pends on the zenith angle and the refractive indexes of both 
sides of the interface. In DMRT-ML, it is obtained by Fres- 
nel’s relationships (e.g. Jin, 1994, p. 59) using the effective 
dielectric constants calculated with Eq. (6). 

At the air-snow interface, the energy coming down from 
the atmosphere 7’“ linos is an input variable given by the user. 
It is assumed isotropic in the current version of DRMT-ML 
but this assumption can be easily changed. The calculation of 
y-atmos f rom meteorological data requires an external model 
such as proposed by Rosenkranz (1998) and Saunders et al. 
(1999). 

At the bottom interface, to account for the diversity of me- 
dia that can be present below snow covers (e.g., soil, ice, lake 
ice, sea ice) and the diversity of electromagnetic modeling 
approaches for the soil, several substratum models are avail- 
able in DMRT-ML and the addition of new models is easy. 
The role of the substratum model is to provide the reflec- 
tion coefficient of the interface. DMRT-ML only takes into 
account the reflection of the coherent wave and neglects the 
diffuse reflections at the bottom interface. It means that the 
roughness of the interface is only partially taken into account. 
In addition to reflecting downwelling energy, the substratum 
is an emitter: the energy emitted and entering into the snow- 
pack is calculated using the temperature of the substratum 
given as an input variable and the emissivity deduced from 
the reflection coefficient and Kirchhoff’s law. 

DMRT-ML proposes 1 1 models that cover soil, ice, lake 
and semi-infinite snowpack. The user selects the type of sub- 
stratum model using an identification number and provides 
the required variables and parameters depending on the type 
of model (Table 1). Since these substratum models have been 
published elsewhere, their evaluation is not addressed in this 
paper. In the case where snow lies over ice, DMRT-ML offers 
two options. If the ice is free of bubbles, isothermal and semi- 
infinite, the best option to use is the “ice” substratum model 
(ID 202, or ID 1 with ice dielectric constant). In any other 
conditions, it is recommended to represent the underlying 
ice using layers with a high density (up to 917 kg m -3 for 
bubble-free ice) and no substratum. The number of layers 
to use depends on the temperature and density profile in the 
ice. The total depth of the ice layers must be large enough to 
avoid “leakage” at the bottom. Special attention is required at 
low frequencies (like L band, 1.4 GHz) since ice absorption 
can be very weak and the layers well below 100 m can sig- 
nificantly contribute. 
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Table 1 . Available models for the reflection coefficient of the bottom interface and required input variables and parameters. SM is the soil 
moisture (volume fraction), s is the complex dielectric constant, a is the surface root mean square height (in meters), f c i ay and /sand are the 
fractions of clay and sand, Porga is the density of dry organic matter (kg m -3 ), T is the temperature (in K) and Q and H are dimensionless 
parameters. P99 stands for Pulliainen et al. (1999), D85 for Dobson et al. (1985), M87 for Matzler and Wegmiiller (1987) and M06 for 
Matzler et al. (2006). 


Interface model 

Material 

Dielectric constant 

ID 

Variables and parameters 

No interface 

None 


0 


Flat surface, Fresnel 

Any 

prescribed 

1 

£ 

Flat surface, Fresnel 

Soil 

from P99 

2 

SM, /day? /sand? Porga? T 

Flat surface, Fresnel 

Soil 

from D85 

3 

SM, T 

Rough surface WM99 

Any 

prescribed 

101 

a,s 

Rough surface WM99 

Soil 

from P99 

102 

cr, SM, /day, ,/sand> Porga, T 

Rough surface WM99 

Soil 

from D85 

103 

a, SM, T 

Flat surface, Fresnel 

Ice 

from M87, M06 

202 

T 

QH model W83 

Any 

prescribed 

301 

a, Q, H, s 

QH model W83 

Soil 

from P99 

302 

° r ? Q? H, SM, /day? /sand? Porga? T 

QH model W83 

Soil 

from D85 

303 

a, Q, H, SM, T 

Flat surface, Fresnel 

Fresh water 

from M87 

402 

T 


2.5 Solution of the radiative transfer equation using the 
DISORT method 

Once the extinction and scattering coefficients of every layer 
and the reflection coefficients at all the interfaces are known, 
the radiative transfer equation is solved using the DISORT 
method (Chandrasekhar, 1960). This method takes into ac- 
count multiple scattering within and between the layers, 
which is an asset with respect to the iterative method (Tsang 
et al., 1985; Jin, 1994; Ishimaru, 1997) for which the number 
of calculated order of scattering is limited. It also computes 
the energy propagation in an unlimited number of directions 
(or “streams”) as opposed to the two-stream (Pulliainen et al., 
1999) or six-stream (Wiesmann and Matzler, 1999) methods 
whose formulations are based on a fixed and small number 
of directions. The drawback of the DISORT method is usu- 
ally the computational cost. However, in the case of passive 
remote sensing, isotropic media, and when the phase func- 
tion has a simple analytical form, the azimuthal dependence 
in Eq. (1) can be analytically integrated. This simplifies the 
equation and reduces the numerical complexity and compu- 
tation cost with respect to other cases like active remote sens- 
ing (Stamnes et al., 1988; Picard et al., 2004) or anisotropic 
media. For snow passive microwave modeling, the assump- 
tion of an isotropic medium is reasonable and the DMRT the- 
ory with small scatterers predicts that the phase function is 
the Rayleigh phase function (Jin, 1994, Eqs. 2-26 to 2-28), 
which only involves the cosines of the azimuth angle. 

In the single layer formulation of DISORT, the integra- 
tion over the zenith angle O' uses a Gaussian quadrature (Jin, 
1994, p. 96), i.e., it is replaced by a discrete sum of inte- 
grand evaluation at optimal angles. The number of angles n 
is defined by the user. This approach cannot be seamlessly 
transposed for multiple layers. The variations of refractive 


index - mainly driven by the snow density profile - cause 
a change of the direction of the streams between the layers 
(Fig. 2). A possible approach (Liang et al., 2008) uses the 
same Gaussian quadrature in every layer (as in the single- 
layer case), which ensures an optimal integration at the ex- 
pense of complex boundary conditions since cubic spline in- 
terpolations are needed to “reconnect” the streams. We use 
a simpler approach in DMRT-ML similar to that proposed 
in Jin (1994, p. 151). The angles are determined by Gaus- 
sian quadrature only in the most refractive layer. In the other 
layers, the angles are determined from Snell’s refraction law, 
which ensures the continuity of the streams between the lay- 
ers. The boundary conditions are simpler at the expense of a 
sub-optimal integration (except in the most refractive layer). 
This issue is easily compensated by increasing the number of 
streams n. 

Another consequence of using the refraction law to deter- 
mine the stream angles is that the number of streams (n, t, 
k = 1 • ■ • L, where L is the number of layers) varies from one 
layer to another. This is caused by the total reflection at large 
zenith angles when a stream propagates from a higher to a 
lower refractive index layer. Such streams in the high refrac- 
tive index layer have no companion in the low index one. 
This also applies to the no streams emerging from the snow 
into the air. Since snow density is much higher than air den- 
sity, no is usually much lower than n. The only consequence 
for the user is that n is not the number of emerging streams as 
it would be in the six-stream or two-stream methods. In prac- 
tice, it is recommended to adjust n to get the wanted number 
of emerging streams no or alternatively to increase n until 
the residual variations of brightness temperature are less than 
the wanted accuracy. Figure 2 illustrates a 4-layer snowpack 
(with snow density values, from top to bottom, of 50, 400, 
200, 320 kgm -1 ). Only upwelling streams are represented. 
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Fig. 2. Upwelling streams for a 4-layer snowpack with varying 
density. The stream directions are calculated with Snell’s refraction 
law using the real part of the effective propagation constant 
(Eq. 6) that mostly depends on snow density. The number of streams 
is 8 in the most dense layer and decreases as the density decreases. 
Only 4 streams emerge in the air above the snowpack. For clarity, 
downwelling streams are not represented. 


In this particular example with extreme variations of density, 
the number of streams in the air is only 4 whereas it is 8 in the 
most refractive layer. The total reflection phenomenon is also 
called internal reflection in Wiesmann and Matzler (1999) 
and causes energy trapping when a layer is surrounded by 
less dense layers. 

Once the integration is replaced by the discrete sum, the 
differential Eq. (1) in z is written for every angle 9: (j — 
±1 ■ ■ ■ ± nk, j >0 and j < 0 for the upward and downward 
directions, respectively) in each layer. 

d 

lij — T H (jij,z) = -K e T B (Hj,z) 

in 

+ £ wyP(/r j , f! / v)r B (M , / - , z) + K a T I (9) 

r 

where jXj — cos (9/), P is found in Jin (1994, Eq. 2-28) and 
I is the unit column vector. In the following, we only present 
the most relevant equations, intermediate calculations are 
given in Jin (1994, Chaps. 4 and 5). 

By introducing the asymmetric and symmetric bright- 
ness temperatures — T B (fij) + Tb( — dj) and 

TBa(Vj) = respectively), Eq. (9) be- 

comes 


l^j — T Ba (fj,j,z) = -K e T Ba (jij,z), 

d " k 

Hj-^T Bs (nj,z) = £ a jji T Ba (At'/ , z) + 2y a T, 


Ajj' = —K e lSjji + -K t wy 


[2(l-^)(l_^,) + ^2] p 


(10a) 

(10b) 

(10c) 


where I is the 2x2 identity matrix and <5 , is 1 when j — j', 
otherwise 0. Note that the W matrix in Jin (1994, Eq. 4-42b) 
is trivial in our case W = —K e owing to the symmetry of the 
Rayleigh phase function. Differentiating Eq. (10a) and com- 
bining it with Eq. (10b) leads to a second-order ordinary sys- 
tem of differential equations. The solutions are of the form 

T Ba(M;,Z) 


2 rii 

=£[* m cxpfo'j,; £)Ty m cxp( o/ m (z T d))J Tg a m 2 TI , (11) 

m 

where d is the layer depth, a m — +«/A m , (m — 1 ...2n , ) 
and A m are the eigenvalues and T\[ a m (Hj) the eigenvec- 
tors of the 2/7, x 2 n,- matrix whose elements are ^kA r ' p , 

11 

(j = p = v,h and j' — p' — v, h). x m and y m 

are 2x2 n, unknown constants to be determined with the 
boundary conditions. In DMRT-ML, the eigenvalue problem 
is solved using LAPACK routines (Anderson et al., 1999). 

Combining Eqs. (10b) and (11), the solutions for the sym- 
metric brightness temperature are 


2 Hj 

T Bs (Hj,z) — £ [x m exp(a m z) (12a) 

m 

-y,n exp (-«,„ (z + d))] T^ s m (nj), 

n k j 

T° Bs ,M = £ —A jj'T^Qij,) (12b) 

j / 


Note that T Bam (dj) and T B% m (dj) are written E and Q, 
respectively, in Jin (1994, p. 102). 

The boundary conditions express the conservation of en- 
ergy at every interface. For each layer k of depth dk, every di- 
rection j and polarization, there are two boundary conditions 
(Jin, 1994, Eqs. 5-10a and 5-10c). The boundary condition 
at the top interface is 


x m,k 


m 

+ [( J ~ r ^fc) r Ba ,m,k^hk)+ 

(l + r Tk) r Bs exp (-a m , k dk)y m ,k 

n k - 1 

m 

\ (M j,k) 1 ( &m,k—ldk—l) x m,k—\ 


r 0 

Bs, m,k— 1 


Ba ,m 


,k-\ 


+ T Bs, 

= (l-r ^(Tk-i-Tk), 


y m ,k - 1 


(13) 
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and the one at the bottom interface is 
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= 1 


(l — r^ffi-TWO, 


(14) 


where r 4 °^ and r^ l, “ onl are the reflection coefficients for 
streams going from the layer k to the upper k — 1 and lower 
k + 1 layers respectively. At the top of the snowpack, Eq.( 13) 
is formally valid if xq and yo are set to zero and 7o is T B tmos . 
At the bottom of the snowpack, Eq. (14) is valid if x^ + \ and 
yL+i are set to 0, the reflection coefficients r^ Uom are cal- 
culated with the substratum model (Table 1) and 7/,+i is the 
substratum temperature. 

The set of boundary conditions forms a linear system of 
N equations ( N — 4 n/.). Since the unknowns x and y in 

layer k are only linked to unknowns in layer k — 1 and k+ 1 , 
the system can be arranged in block-diagonal structure and 
can be solved using the efficient banded matrix solver in LA- 
PACK. The brightness temperature emerging from the snow- 
pack is then calculated with 





[^Ba,m, lOf/,o) + ^Bs.m.lt/T/.o)] x m, 1 
+ [^Ba.m.l ~ ^Bs,m,l(/T/d) 


exp ( (Xm . i d [ ) y m , i ) • 


(15) 


The value from this last equation is returned to the user as 
the simulated top-of-snowpack brightness temperature. 


2.6 Post-computation: emissivity and reflectivity 

The emissivity, i.e., the coefficient measuring the departure 
from a black body, is often used to present modeling results 
in passive microwave studies instead of brightness temper- 
ature (whose value is related to the snow physical tempera- 
ture). Only when the medium is strictly isothermal, that is, 
when the physical temperature of snow and the underlying 


medium is uniform and equal to T, it is possible to compute 
the emissivity using a single simulation and the following 
equation: 



For a nonconstant temperature profile, which is the rule 
for any natural snowpack, the definition of emissivity is not 
trivial. To mimic the in-equilibrium case and the Kirchhoff 
law, the emissivity can be defined as one minus the reflec- 
tivity of the medium (Peake, 1959). The calculation requires 
two simulations with slightly different atmospheric bright- 
ness temperatures (:rg tmos and 7^ tmos + Arg tmos ). Assuming 
T b and T' B are the results of these simulations, the reflectiv- 
ity and emissivity are given by 


e—l — r = 1 — 


T'b-T b 

a yatmos * 


(17) 


In practice, using 7’j“ mos = 0 and Arg tmos = 1 K is recom- 
mended. 


3 Results 

This section presents the sensitivity of DMRT-ML to the 
most important snow properties required as inputs. It also 
discusses the limitations of the model and provides the range 
of validity of the input variables whenever possible. 

3.1 Sensitivity to ice dielectric constant 

Figure 3a shows brightness temperatures at 1 8 GHz as a func- 
tion of the zenith angle calculated with DMRT-ML and cal- 
culated by Kong et al. (1979) together with observations re- 
ported by Tsang and Kong (2001, Fig. 7.7.2). The medium 
is considered semi-infinite, the grain size was chosen to be 
1.75 mm and the ice dielectric constant £j ce = 3.2 + / 0 . 0 1 6 . 
The result of the DMRT-ML simulation with prescribed di- 
electric constant (solid line) closely matches the original cal- 
culation (dotted line) by Tsang and Kong (2001) up to in- 
cidence angles of 65°. This provides a technical validation 
of our implementation in the single layer case. However, 
the imaginary part of the ice dielectric constant used by 
Kong et al. (1979) was an order of magnitude higher than 
the one obtained with the more recent Eq. (2) from Matzler 
and Wegmiiller (1987). With the latter parameterization, the 
brightness temperatures simulated using the same grain size 
of 1.75 mm are much lower (dash line in Fig. 3b), leading to a 
strong disagreement with the observations. Modeling results 
and observations can be re-conciliated by using a smaller ra- 
dius of 0.83 mm (solid line in Fig. 3b). This new simulation 
yields results very close to those of the original simulation 
and observations. 

Even if Eq. (2) is likely closer to reality than earlier for- 
mulas, the parameters of ice dielectric constants are still not 
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Fig. 3. Comparison of DMRT-ML simulations (solid line) at 18 GHz with previous calculations (Tsang and Kong, 2001) (dotted curve) 
and with experimental observations (filled dots). The medium is semi-infinite, with a density of 350kgm -3 , a temperature of 272 K, a 
grain radius of 1.75 mm and an ice dielectric constant of 3.2 + ;0.016. Simulations with a more realistic dielectric constant (Matzler and 
Wegmiiller, 1987) are shown on the right panel with the original radius of 1 .75 mm (dashed curve) and refitted radius of 0.83 mm (solid line). 
The gray bars represent variations of the dielectric constant imaginary part of ±20 % around in the later case. 


well constrained by measurements especially at low frequen- 
cies (Jiang and Wu, 2004) and the uncertainty is unknown. 
The consequences on the predicted brightness temperature 
can be significant as illustrated in Fig. 3b (gray area) obtained 
by assuming an error of ±20 % of the dielectric loss (imagi- 
nary part of the dielectric constant). At a 53° incidence angle 
for instance, the error is 1 1 and 12 K at vertical and horizon- 
tal polarizations respectively. Such an error is significant and 
must be taken into account in the interpretation of simula- 
tions by thermal microwave emission models. 

3.2 Sensitivity to the grain size 

The significant sensitivity of microwave thermal radiation to 
snow grain size is widely recognized (Zwally, 1977). It stems 
from the fact that snow grains are usually smaller than the 
wavelength in the microwave range and their scattering co- 
efficient (Eq. 8c) increases as the cubic power of the sphere 
radius a. Figure 4 presents the variation of the scattering ef- 
ficiency (2 S = 4cik s /3 /) as a function of sphere radius. The 
calculations with DMRT-ML (for various densities) are plot- 
ted as solid lines. 

In contrast, the absorption coefficient increases linearly 
with the size. It results that the emissivity and the bright- 
ness temperature of a semi-infinite medium strongly decrease 
with size (Fig. 5). The difference between vertical and hori- 
zontal polarizations also slightly decreases as scattering in- 
creases. 

An important consequence of this cubic dependence is that 
a collection of spheres with different sizes is not equivalent to 
a collection of identical spheres with the averaged size. The 
contribution of the largest spheres of the collection to scatter- 
ing is comparably greater than the contribution of the small- 


est. However, Jin (1994) shows that, under the assumption 
of small grains, any collection can be represented by equal- 
radius spheres with an equivalent grain size a c . This size is 
always larger than the average of the distribution but depends 
on the distribution shape as well as snow density (Jin, 1994, 
Eq. 3—42 and Figs. 3-11, 3-12). In DMRT-ML, we imple- 
mented the calculation for a Rayleigh distribution of size and 
reached the same conclusion as Jin ( 1 994). Unfortunately, the 
results are highly dependent on the choice of the distribution 
and especially on the slope of the upper tail of the distribu- 
tion. For instance, a log-normal distribution - relevant for 
snow (Flanner and Zender, 2006) - features many very large 
grains for a reasonable mean grain size, which leads to the 
divergence of the DMRT-ML calculation and breaks the as- 
sumptions of small scatterers. In practice, the choice of the 
distribution is difficult and is related to the more general is- 
sue of the representation of snow by a collection of spherical 
grains. 

Figure 4 also illustrates the empirical correction for large 
particles proposed in DMRT-ML (Sect. 2.3). The corrected 
scattering efficiency (circles) is limited to a maximum value 
of 2, which corresponds to the theoretical asymptotic value 
for very large particles (Grody, 2008). This correction elim- 
inates the unrealistic divergence of the scattering efficiency 
for large grain sizes (solid line). Nevertheless, the quality 
of this correction is difficult to evaluate for dense media. 
For a sparse medium however, the scatterers are independent 
and the scattering efficiency is obtained by Mie’s calculation 
(Warren, 1982). Figure 4 (dashed blue curve) shows that the 
Mie scattering efficiency at 89 GHz diverges from DMRT- 
ML with density tending to 0 kg m -3 (i.e., independent scat- 
terers) for grain radii larger than 0.75 mm (i.e., a A./5 
and Q s of 2-3). This is in agreement with our correction. 


www.geosci-model-dev.net/6/1061/2013/ 


Geosci. Model Dev., 6, 1061-1078, 2013 


1070 


G. Picard et al.: Simulation of the microwave emission of multi-layered snowpacks 



Fig. 4. Scattering efficiency (Q s = 4aK s /3f) at 89 GHz as a func- 
tion of the grain size for different approximations and densities: 
Independent scattering approximation (applies for very small den- 
sity only) using DMRT-ML with a density approaching 0 (solid 
blue line), or full Mie calculation (dashed blue line); DMRT-ML 
QCA-CP nonsticky with a density of 200 kg m -3 (solid red line) 
and 300 kg m -3 (solid black line). Dotted curves show calculations 
based on the DMRT-ML theory with a limitation of the scattering 
efficiency at Q s = 2, for two different densities. 



Grain radius (mm) 

Fig. 5. Brightness temperature at 89 GHz, 53° incidence angle, ho- 
rizontal and vertical polarizations (dashed and solid line, respec- 
tively) as a function of the grain size with a density of 200 (red) and 
300 kg m -0 (black). 

However, Mie efficiency reaches a maximum value of nearly 
5 and remains largely above 2 in the range of realistic grain 
sizes. This is not reproduced by our correction. In fact, the 
convergence towards 2 is only observed at much larger grain 
sizes. This result suggests that the empirical correction would 
be more accurate if the efficiency limit were increased up to 
a value around 4. However, this result was derived in partic- 
ular conditions (sparse medium, 89 GHz) and its generality 
is tin knovvn. We therefore recommend to use the correction 
with caution and only when a very limited number of layers 
in the snowpack have large grains. It is worth noting that the 
condition Q s < 2 is valid for most types of snow at frequen- 
cies lower or equal to 89 GHz as shown in Fig. 6. Harlow 
and Essery (2012, Fig. 1 1) show using Mie-QCA with sticky 
spheres (r = 0.2) that this assumption is valid up to about 
60 GHz for 1 mm-radius particles and to 200 GHz for par- 
ticles of 0.3 mm radius and smaller. This confirms that the 



Fig. 6. Range of grain sizes for which the DMRT-ML QCA-CP non- 
sticky is reasonably valid as a function of the frequency for snow at 
260 K and 200 (red) or 300kgm“° (black). The upper limit is the 
grain size for which the scattering efficiency reaches a value of 2. 


general criteria a ;$ X/5 is adequate to evaluate the validity 
of QCA-CP for small scatterers. 

However, even if the condition Q s < 2 is true, we note that 
the brightness temperatures become unrealistically low be- 
fore Q s reaches a value of 2 (which occurs for a > 1 mm in 
Fig. 5). For example, the lowest brightness temperature ever 
observed at 89 GHz in Antarctica by AMSR-E is 1 17 K. 

3.3 Influence of snow density 

Snow density is involved in many components of the model. 
It drives (i) the proximity of the grains, thus the scatter- 
ing coefficient of the medium in relation with the grain size 
(Eqs. 5, 6); (ii) the mass of ice, thus the absorption coeffi- 
cient; (iii) the real part of the refractive index, which deter- 
mines the stream angles, and the transmission and reflection 
coefficients of every interface (West et al., 1996). 

The first two effects are illustrated in Fig. 7 at 37 GHz and 
for a grain size of 0.3 mm. The absorption and scattering co- 
efficients are calculated assuming a medium of “ice spheres 
in air” for density less than half of the pure ice density (i.e., 
fractional volume of 50%) and “air spheres in ice” other- 
wise. The discontinuity observed in Fig. 7 comes from the 
fact that both representations are not strictly equivalent even 
at a fractional volume of 50 %. 

The absorption coefficient increases almost linearly with 
the density because the ice is the only absorber. In contrast, 
the case of the scattering coefficient is more complex. For 
very low density, the medium is sparse and the scattering 
coefficient calculated with the DMRT-ML theory increases 
closely to the one calculated with the independent scatter- 
ers assumption (Fig. 7, dashed line). However, for densities 
larger than 100 kg m -3 , the latter becomes invalid, because 
the scatterers are too close to be considered independent - 
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they are in fact in the shadow of each other, which weak- 
ens their scattering efficiency (Liang et al., 2006). This effect 
is nicely captured by the DMRT-ML theory (Tsang et al., 
2000b, Fig. 10.4.5), which predicts that the rate of variation 
of the scattering coefficient decreases with density. The co- 
efficient reaches a maximum at a density around 150 kg m -3 
and decreases for higher densities. Although this general be- 
havior is expected, the accuracy of the DMRT-ML theory 
is not well known at such high densities. Recent work us- 
ing exact numerical calculation has shown that the theory 
DMRT-ML QCA starts to deviate from fractional volume 
around 30% (Liang et al., 2006), i.e., a density of 275 kg 
m -3 (validity range represented by circles in Fig. 7). How- 
ever, the generality of this result for smaller grains, moderate 
stickiness or under the QCA-CP assumption is unexplored. 
If we assume that the value of 30% is correct and applies 
also for “air spheres in ice”, the theory would be valid in 
the range 640-917 kg m -3 (represented by squares in Fig. 7). 
Unfortunately, snow density falls in the intermediate range 
(275-640 kg m -3 ) in many conditions. To deal with this is- 
sue in practice, a pragmatic option is to consider the devia- 
tion above 30% fractional volume is moderate with respect 
to other errors (such as grain size measurements) and to apply 
DMRT-ML QCA-CP using the most adequate medium rep- 
resentation as a function of the density as in Fig. 7. A second 
option is to interpolate the scattering and absorption coeffi- 
cients using polynomials fitted with anchor points taken in 
both domains where the theory is valid. This option called 
“bridging” (Dierking et al., 2012) is appealing because it 
yields continuous relationships as a function of the density. 
Flowever, the sensitivity to the choice of the polynomial or- 
der and the anchor points has to be evaluated. Therefore, the 
“bridging” option is not implemented in DMRT-ML yet. 

3.4 Influence of the stickiness 

Figure 8 shows the scattering coefficient at 37 GFIz as a func- 
tion of density for different values of the stickiness parameter 
r and grain radius a for a given temperature. It shows that the 
scattering coefficient increases as the stickiness parameter r 
decreases (from blue to black to green curves). The lower 
values of stickiness correspond indeed to stronger attractions 
between spheres and a more pronounced clustering effect. 
Clusters are “seen” by the microwaves as larger objects than 
the particles they are made of. It means that they scatter more 
than if the particles were randomly positioned (i.e. nonsticky 
case). However, the stickiness parameter is not just a scaling 
factor of the grain size. In fact, a cluster of small particles is 
not equivalent to a large particle as illustrated in Fig. 8: the 
sensitivity to the density is different between a cluster (green 
curve) and a large particle (red curve). The stickiness tends 
to shift the maximum of the scattering coefficient toward a 
larger density. 

In practice, choosing a realistic value of stickiness to rep- 
resent natural snow is difficult. There is currently no means 



Fig. 7. Scattering (blue and plain symbols) and absorption (red and 
hollow symbols) coefficients at 37 GHz as a function of the density. 
The temperature is 260 K and the grain radius is 0.3 mm. The model 
“ice spheres in air” is used for densities less than half the pure ice 
density (458.5 kg m -3 ). For higher densities, the model “air spheres 
in ice” is used. Circles and squares show the domain of validity of 
the DMRT-ML theory for each model. The scattering coefficient 
under the assumption of independent “ice spheres in air” is shown 
as a blue dash curve. 



Fig. 8. Scattering coefficient at 37 GHz as a function of the density 
for different stickiness parameters r and grain radius a. The tem- 
perature is 260 K. 

to estimate this value either from field measurements, 3-D 
images of snow micro-structure or snow evolution model 
outputs. As for choosing a grain size distribution, the core 
of the problem is the representation of snow by spheres. 
Tsang et al. (2008) suggest to use r = 0.1 because it yields 
a frequency-dependence in better agreement with measure- 
ments. With such a low value, the clustering effect is in 
fact very strong and the size of the cluster approaches the 
wavelength (long-range effect), which explains the change 
of the frequency-dependence. Matzler (1998) suggests to use 
a higher value, r = 0.2, based on a comparison between the 
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density dependence in DMRT-ML with various degrees of 
stickiness and in the improved Born approximation (Matzler, 
1998). 

The implementation of the stickiness in DMRT-ML is lim- 
ited to the “short range” version, i.e., both the grains and the 
cluster must be small with respect to the wavelength. In prac- 
tice, r should be larger than its theoretical minimum value 
of (2 - V 2)/6 = 0.098 (Tsang et al., 2000b, p. 427). After 
Tsang et al. (2000b, Fig. 8.4.3) and our own calculation (not 
shown), the “short range” calculation is valid for r values 
larger than 0.25 and grain sizes, which respects the small 
scatterer assumption. For lower r, the validity depends on 
the grain size. 

3.5 Influence of the liquid water content 

Figure 9 shows the relationship between the brightness tem- 
perature at 19 GHz and horizontal polarization as a func- 
tion of the total column liquid water content. The snowpack 
is considered homogeneous, except that the liquid water is 
concentrated in the top 10 cm. This calculation confirms the 
strong influence of the liquid water on brightness tempera- 
ture, which is exploited to detect snowmelt events from pas- 
sive microwave observations (e.g. Picard and Fily, 2006). It 
also shows that brightness temperature reaches a nearly con- 
stant value from about 0.5 kgm -2 of liquid water. This fea- 
ture explains why the retrieval of the amount of liquid wa- 
ter from passive microwave observation is probably impos- 
sible. The threshold value, 0.5 kgm -2 , is close to value ob- 
tained with the MEMLS model (Tedesco et al., 2007) and 
that obtained indirectly by comparing observations and out- 
puts of the RACMO regional snow and atmosphere model 
(Kuipers Munneke et al., 2012). 

3.6 Layered snowpack 

It is well known that the natural snowpack is composed of 
relatively homogeneous layers formed by episodes of accu- 
mulation and subject to metamorphisms. Hence, the density 
and grain size (and stickiness) can be different in every layer. 
Since the electromagnetic properties (i.e., the scattering and 
absorption coefficients and the effective constant of propaga- 
tion) have a nonlinear dependence to the snow properties, the 
single-layer “average” snowpack is not electromagnetically 
equivalent to the multi-layer snowpack. In addition, the tem- 
perature is rarely uniform within the snowpack, especially 
close to the surface where the strong temperature gradients 
are caused by the daily variations of solar energy (Brandt 
and Warren, 1993). 

Several effects can result from the layered nature of the 
snowpack. In general, accounting for the multi-layered na- 
ture is particularly important for 

- the difference between the horizontal and vertical po- 
larizations. The difference is particularly sensitive to re- 
flection at the air-snow interface, which is driven by the 



Fig. 9. Brightness temperature at 19 GHz and horizontal polariza- 
tion as a function of the total amount of liquid water in the snow- 
pack. The temperature is 273 K, density is 300 kgm -3 and there 
exist two grain size values (0.5 and 1 mm). The liquid water is con- 
centrated in the top first 10 cm of the snowpack. 

density in the upper layer and by the reflection at the 
internal interface, which depend on the contrast of den- 
sity between successive layers. There are several exper- 
imental and theoretical evidences of this effect (Matzler 
et al., 1984; Liang et al., 2008; Champollion et al., 
2013). Smoothing the profile of density in simulation 
results directly in a decrease of the difference between 
the polarizations. 

- the spectrum of emissivity. It is particularly sensitive to 
the scattering and absorption profiles because the pen- 
etration depth highly depends on the frequency. For in- 
stance, Bmcker et al. (2010) showed that the spectra ob- 
served in Antarctica cannot be explained with a homo- 
geneous snowpack, and that the increase of grain size 
with depth is a necessary condition to explain the ob- 
served spectra. In some extreme cases called “anoma- 
lous spectra” by Rosenfeld and Grody (2000), the ob- 
served emissivity increases with frequency although the 
emissivity of a homogeneous snowpack decreases with 
frequency due to a stronger increase of scattering rel- 
ative to that of the absorption. To illustrate the effect 
of the resolution of the snow parameter profiles, the 
calculation presented in Brucker et al. (2011a) using 
the original 2.5 cm high-resolution profiles of grain size 
and density measured at Dome C in Antarctica are pre- 
sented in Fig. 10 along with calculations using lower- 
resolution profiles. The simulations are performed at 
19 and 37 GHz with a homogeneous temperature of 
219 K (the annual mean at Dome C). Only the resolution 
in the upper 3 m for which measurements were avail- 
able is varied. Figure 10 shows that coarser-resolution 
results in higher brightness temperatures (up to about 
8 %), except for the lowest resolution (3 m-thick layer, 
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corresponding to the homogeneous snowpack). This 
complex dependence might be explained by the fact that 
the modeled brightness temperature is more sensitive to 
the resolution of the density profile than that of the grain 
size (simulations not shown) and that the density depen- 
dence of the electromagnetic properties is not mono- 
tonic as shown in Fig. 7. These results are specific to 
Dome C but they emphasize the importance of the reso- 
lution of the input parameters. 

- time series of brightness temperature. Even if the grain 
size and density are assumed homogeneous the temper- 
ature profile is rarely uniform nor constant over time in 
snow and its variations can even be significant in the 
layer identified in the field based on the homogeneity of 
the grains and compactness. The variations of tempera- 
ture near the surface are in general more rapid than the 
change of grain size and density due to metamorphism, 
and thus have the most important contribution to the 
short-term brightness temperature variations. Simula- 
tions at Dome C in Brucker et al. (20 11a) show that even 
with grain size and density profiles taken constant over 
a few years (metamorphism is very slow at Dome C due 
to low temperature), the variations of brightness temper- 
ature are well reproduced (^ 2 K) using measured time 
series of temperature profiles. Reducing the resolution 
of the temperature profile would result in a smoother 
time series. 



Fig. 10. Brightness temperature at 18.7 and 36.5 GHz at Dome C, 
Antarctica, simulated as a function of the resolution of the pro- 
files of grain size and density. The original profiles are presented 
in Brucker et al. (201 la) and are composed of two parts: the upper 
3 m were measured with a vertical resolution of 2.5 cm (simulation 
marked by a symbol) and the lower part (3 to 100 m depth) is a de- 
terministic function of the depth. The profiles at lower resolutions 
are generated by merging successive layers in the upper part only 
to emphasize the influence of measured profile resolution, the lower 
part remaining unchanged. The x axis reports the thickness of the 
layers, from 2.5 cm (original profile) to 3 m (i.e. one single layer for 
the upper part). 


4 Validation of DMRT-ML with external data 

The comparison between measured brightness temperatures 
and results of DMRT-ML simulations using measured in- 
puts was addressed in several studies: for a typical dry semi- 
infinite snowpack in Antarctica (Brucker et al., 2011a), for 
Arctic and sub-Arctic seasonal snowpacks (Roy et al., 2013) 
and for snow overlying ice as found in the ablation zone of 
the ice sheets (Dupont et al., 2012). In the three cases, it was 
necessary to estimate some parameters by optimization with 
respect to the measured brightness temperature. It is indeed 
difficult to measure all the input variables and parameters 
and the brightness temperatures with the same representa- 
tiveness. In addition, some quantities - like the snow grain 
size and the soil properties - are notoriously difficult to de- 
termine in the field. The representation of snow grain in the 
DMRT-ML theory by spherical particles is also a conceptual 
difficulty. 

Hence, the results of the comparisons and the errors esti- 
mated in these studies depends on the methodology and are 
meaningless out of the context of each study. Nevertheless, 
these studies converge on two facts. First, the grain size de- 
rived from specific surface area measurements, i.e., optical 
radius, needs to be increased by a factor between 2.8 and 
3.5 to be suitable as input of DMRT-ML (see discussion in 


Roy et al., 2013). There is no evidence that similar adjust- 
ments would be required for the density or temperature mea- 
surements. Second, the model predicts reliable dependence 
between horizontal and vertical polarizations near the Brew- 
ster angle (50-55°). Figures 11 and 12 illustrate the latter 
point for the polarizations at 37 and 19 GHz respectively. 
They show the brightness temperature at horizontal polariza- 
tion versus vertical polarization for all the snow pits or pixels 
analyzed in Brucker et al. (2011a), Roy et al. (2013), and 
Dupont et al. (2012). Observations (in red) were acquired 
with a ground-based radiometer except at Dome C (stars), 
where measurements were recorded by the advanced mi- 
crowave scanning radiometer for EOS (AMSR-E). DMRT- 
ML predictions are in blue. A large variety of environments 
are represented: Antarctica, ice-sheet ablation area, Arctic 
tundra, Arctic windy tundra, Arctic fen, and grassland. Each 
gray line links the observation and simulation result from the 
same snow pit or pixel. The length of each line illustrates the 
discrepancy between the observation and the simulation re- 
sult of the order of 2-13 K but, as stated before, this is not an 
absolute error since it depends on the calibration procedure 
that differs between the studies. The relationship between the 
polarizations can be safely interpreted as it is almost inde- 
pendent of the calibration: Figs. 11 and 12 demonstrate that 
the measured brightness temperatures at both polarizations 
are highly correlated over a large range of about 120 K at 
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Fig. 11. Brightness temperatures at 55° incidence angle and 
37 GHz, at vertical polarization versus horizontal polarization for 
a variety of sites: Dome C, Antarctica (stars), Arctic tundra (trian- 
gles up), Arctic windy tundra (triangles down), Arctic fen (squares), 
grassland (diamonds), ice-sheet ablation or percolation areas (cir- 
cles). Observations are in red and DMRT-ML predictions in blue. 
The gray lines link the observation and prediction of the same site. 




Fig. 12. Same as Fig. 11 for a frequency of 19 GHz. 


37 GHz and 70 K at 19 GHz. Part of this correlation stems 
from the linear relationship between the brightness tempera- 
ture and the snow physical temperature. However the range 
of physical temperature - typically from —55 to 0°C - is 
less than the brightness temperature ranges. It means that the 
emissivities at both polarizations are correlated. At 19 GHz, 
the correlation is lower than at 37 GHz. A probable explana- 


tion is the larger contribution of the substratum at the lower 
frequency due to the longer penetration depth. The important 
point is that the model nicely reproduces this general corre- 
lation. Even where the model and the observations disagree 
(i.e. long gray lines), the correlation between polarizations 
remains (i.e. the gray lines follow the general trend of the 
points). 


5 Conclusions 

The DMRT-ML is a physically based model used to compute 
brightness temperature at any frequency in the microwave 
range and at horizontal and vertical polarizations from input 
variables describing multi-layered snowpack and its environ- 
ment. These variables and parameters include the profiles of 
snow temperature, density, grain size, stickiness, and liquid 
water content, the characteristics of the substratum (e.g. soil 
moisture, texture and temperature in the case of a soil sub- 
stratum), and the downwelling atmospheric brightness tem- 
perature. 

The paper presents the sensitivity of the microwave emis- 
sion to the most important input variables and parameters 
and makes recommendations on the validity ranges of these 
variables and parameters, either constrained by the underly- 
ing DMRT-ML theory or by the specific DMRT-ML imple- 
mentation. The validation of DMRT-ML with external in situ 
measurements is detailed in several studies (Brucker et al., 
2011a; Roy et al., 2013; Dupont et al., 2012) for various en- 
vironments. The error found between predicted and observed 
brightness temperatures ranged between 2 and 13 K, which 
gives the magnitude of accessible errors but which depend 
on the methodology used for the comparison. In particular, 
the choice of the relationship to relate the measured grain 
size to the grain size metric relevant to the DMRT-ML the- 
ory is critical and no ideal solution exists yet. However, these 
studies add up to many others that have contributed to vali- 
date the DMRT-ML theory (Macelloni et al., 2001; Tsang 
and Kong, 2001; Tedesco et al., 2004; Grody, 2008). Cur- 
rently, the most problematic point is probably the limited ac- 
curacy of the DMRT-ML theory for intermediate density val- 
ues (about 300-500 kg m -3 ) that are commonly observed in 
natural snow and firn. An empirical correction of this prob- 
lem has been recently proposed (Dierking et al., 2012). Even 
though it could be accurate enough, it has not the theoretical 
grounds that characterize the DMRT-ML theory and makes 
one of its merits. Further improvements are needed. 

The main characteristics of the DMRT-ML implemen- 
tation are the availability as an open source software, the 
efficiency of the computation due to the use of Fortran90 
and LAPACK library, the wide range of cryospheric environ- 
ments that can be modeled without any change of the source 
code (dry and wet snowpack, seasonal snow over soil, sea- 
sonal snow over bubbly ice, frozen lakes, perennial snow, 
etc.). In addition, modeling sea ice will be possible in the near 


Geosci. Model Dev., 6, 1061-1078, 2013 


www.geosci-model-dev.net/6/1061/2013/ 


G. Picard et al.: Simulation of the microwave emission of multi-layered snowpacks 


1075 


future with the implementation of the dielectric constants 
of salted snow and water. These particularities are strong 
assets for the coupling of DMRT-ML with land surface 
models (Vionnet et ah, 2012) and for the integration as an 
observation operator into data assimilation schemes or com- 
putationally intensive inverse methods. DMRT-ML is avail- 
able for download at http://lgge.osug.fr/~picard/dmrtml/. 
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